MODULE FFT
  USE Precision
  IMPLICIT NONE
  PRIVATE
  PUBLIC :: FFTIntegrate,DFTCorrection
CONTAINS
  SUBROUTINE FFTIntegrate(func,a,b,w,sinint,cosint,params)
    USE Parameters
    USE PolynomialInterpolation
    INTERFACE
       FUNCTION func(x,params)
         USE Precision
         REAL(d) :: func
         REAL(d),INTENT(IN) :: x
         REAL(d),INTENT(IN),OPTIONAL :: params(:)
       END FUNCTION func
    END INTERFACE
    REAL(d),INTENT(IN) :: a,b,w
    REAL(d),INTENT(OUT) :: cosint,sinint
    REAL(d),INTENT(IN),OPTIONAL :: params(:)
    INTEGER,PARAMETER :: M=8192 &
                        ,NDFT=4*M &
                        ,MPOL=6
    INTEGER :: j,nn
    REAL(d) :: cdft,cerr,corfac,corim,corre,en,sine,cosine,sdft,serr &
              ,cpol(MPOL),spol(MPOL),xpol(MPOL)
    LOGICAL,SAVE :: init=.TRUE.
    REAL(d),SAVE :: aold=-1.0e30_d &
                   ,bold=-1.0e30_d &
                   ,delta,data(NDFT),endpts(8)
    initialize: IF(init.or.a.ne.aold.or.b.ne.bold) then
       init=.FALSE.
       aold=a
       bold=b
       delta=(b-a)/M
       DO j=1,M+1
          with_params: IF(PRESENT(params)) THEN
             data(j)=func(a+(j-1)*delta,params)
          ELSE with_params
             data(j)=func(a+(j-1)*delta)
          END IF with_params
       END DO
       DO j=M+2,NDFT
          data(j)=0.0_d
       END DO
       DO j=1,4
          endpts(j)=data(j)
          endpts(j+4)=data(M-3+j)
       END DO
       CALL RealFT(data,.FALSE.)
       data(2)=0.0_d
    ENDIF initialize
    en=w*delta*NDFT/(2.0_d*PI)+1.0_d
    nn=min(max(int(en-0.5_d*MPOL+1.0_d)),NDFT/2-MPOL+1)
    DO j=1,MPOL
       cpol(j)=data(2*nn-1)
       spol(j)=data(2*nn)
       xpol(j)=nn
       nn=nn+1
    END DO
    CALL Interpolate(xpol,cpol,en,cdft,cerr)
    CALL Interpolate(xpol,spol,en,sdft,serr)
    CALL DFTCorrection(w,delta,a,b,endpts,corre,corim,corfac)
    cdft=cdft*corfac+corre
    sdft=sdft*corfac+corim
    cosine=COS(w*a)
    sine=SIN(w*a)
    cosint=delta*(cosine*cdft-sine*sdft)
    sinint=delta*(sine*cdft+cosine*sdft)
  END SUBROUTINE FFTIntegrate
!C  (C) Copr. 1986-92 Numerical Recipes Software D041&0(9p#3.

  SUBROUTINE DFTCorrection(w,delta,a,b,endpts,corre,corim,corfac)
    USE Parameters
    REAL(d),INTENT(IN) :: w,delta,a,b
    REAL(d),INTENT(OUT) :: endpts(8),corre,corim,corfac
    REAL(d),PARAMETER :: CROSSOVER=1.0_d ! (numerical_recipes: 5e-2)
    REAL(d) :: a0i,a0r,a1i,a1r,a2i,a2r,a3i,a3r,arg,cosine,cl,cr,sine,sl,sr &
              ,t2,t4,t6,cth,ctth,spth2,sth,sth4i,stth,th,th2,th4,tmth2,tth4i
    th=w*delta
    IF(a.ge.b.or.th.lt.0.0_d.or.th.gt.PI) THEN
       WRITE(0,*)'ERROR: bad arguments to DFTCorrection'
       WRITE(0,*)'omega:',w,'Delta:',delta,'theta:',th
       STOP
    ENDIF
    IF(ABS(th)<CROSSOVER) THEN ! Use series expansion
       t2=th*th
       t4=t2*t2
       t6=t4*t2
       corfac=1.0_d-(11.0_d/720.0_d)*t4+(23.0_d/15120.0_d)*t6
       a0r=(-2.0_d/3.0_d)+t2/45.0_d+(103.0_d/15120.0_d)*t4 &
            -(169.0_d/226800.0_d)*t6
       a1r=(7.0_d/24.0_d)-(7.0_d/180.0_d)*t2+(5.0_d/3456.0_d)*t4 &
            -(7.0_d/259200.0_d)*t6
       a2r=(-1.0_d/6.0_d)+t2/45.0_d-(5.0_d/6048.0_d)*t4+t6/64800.0_d
       a3r=(1.0_d/24.0_d)-t2/180.0_d+(5.0_d/24192.0_d)*t4-t6/259200.0_d
       a0i=th*(2.0_d/45.0_d+(2.0_d/105.0_d)*t2-(8.0_d/2835.0_d)*t4 &
            +(86.0_d/467775.0_d)*t6)
       a1i=th*(7.0_d/72.0_d-t2/168.0_d+(11.0_d/72576.0_d)*t4 &
            -(13.0_d/5987520.0_d)*t6)
       a2i=th*(-7.0_d/90.0_d+t2/210.0_d-(11.0_d/90720.0_d)*t4 &
            +(13.0_d/7484400.0_d)*t6)
       a3i=th*(7.0_d/360.0_d-t2/840.0_d+(11.0_d/362880.0_d)*t4 &
            -(13.0_d/29937600.0_d)*t6)
    else ! use trigonometric formulas
       cth=cos(th)
       sth=sin(th)
       ctth=cth**2-sth**2
       stth=2.0_d*sth*cth
       th2=th*th
       th4=th2*th2
       tmth2=3.0_d-th2
       spth2=6.0_d+th2
       sth4i=1.0_d/(6.0_d*th4)
       tth4i=2.0_d*sth4i
       corfac=tth4i*spth2*(3.0_d-4.0_d*cth+ctth)
       a0r=sth4i*(-42.0_d+5.0_d*th2+spth2*(8.0_d*cth-ctth))
       a0i=sth4i*(th*(-12.0_d+6.0_d*th2)+spth2*stth)
       a1r=sth4i*(14.0_d*tmth2-7.0_d*spth2*cth)
       a1i=sth4i*(30.0_d*th-5.0_d*spth2*sth)
       a2r=tth4i*(-4.0_d*tmth2+2.0_d*spth2*cth)
       a2i=tth4i*(-12.0_d*th+2.0_d*spth2*sth)
       a3r=sth4i*(2.0_d*tmth2-spth2*cth)
       a3i=sth4i*(6.0_d*th-spth2*sth)
    endif
    cl=a0r*endpts(1)+a1r*endpts(2)+a2r*endpts(3)+a3r*endpts(4)
    sl=a0i*endpts(1)+a1i*endpts(2)+a2i*endpts(3)+a3i*endpts(4)
    cr=a0r*endpts(8)+a1r*endpts(7)+a2r*endpts(6)+a3r*endpts(5)
    sr=-a0i*endpts(8)-a1i*endpts(7)-a2i*endpts(6)-a3i*endpts(5)
    arg=w*(b-a)
    cosine=cos(arg)
    sine=sin(arg)
    corre=cl+cosine*cr-sine*sr
    corim=sl+sine*cr+cosine*sr
  END SUBROUTINE DFTCorrection
!C  (C) Copr. 1986-92 Numerical Recipes Software D041&0(9p#3.

  SUBROUTINE RealFT(data,flagInverse)
    USE Parameters
    REAL(d),INTENT(INOUT) :: data(:)
    LOGICAL,INTENT(IN),OPTIONAL :: flagInverse ! .TRUE. for inverse FT
    INTEGER :: n
    LOGICAL :: doInverse=.FALSE.
    INTEGER :: i,i1,i2,i3,i4,n2p3
!DEREK    REAL(d) :: c1,c2,h1i,h1r,h2i,h2r,wis,wrs
    REAL(d) :: c1,c2,h1i,h1r,h2i,h2r
    REAL(d) :: theta,wi,wpi,wpr,wr,wtemp
    n=SIZE(data)
    IF(PRESENT(flagInverse)) doInverse=flagInverse
    theta=PI/REAL(n/2,d)
    c1=0.5_d
    if (doInverse) then
       c2=0.5_d
       theta=-theta
    else
       c2=-0.5_d
       call DFT(data,.FALSE.)
    endif
    wpr=-2.0_d*sin(0.5_d*theta)**2
    wpi=sin(theta)
    wr=1.0_d+wpr
    wi=wpi
    n2p3=n+3
    DO i=2,n/4
       i1=2*i-1
       i2=i1+1
       i3=n2p3-i2
       i4=i3+1
!DEREK
!       wrs=sngl(wr)
!       wis=sngl(wi)
       h1r=c1*(data(i1)+data(i3))
       h1i=c1*(data(i2)-data(i4))
       h2r=-c2*(data(i2)+data(i4))
       h2i=c2*(data(i1)-data(i3))
!       data(i1)=h1r+wrs*h2r-wis*h2i
!       data(i2)=h1i+wrs*h2i+wis*h2r
!       data(i3)=h1r-wrs*h2r+wis*h2i
!       data(i4)=-h1i+wrs*h2i+wis*h2r
       data(i1)=h1r+wr*h2r-wi*h2i
       data(i2)=h1i+wr*h2i+wi*h2r
       data(i3)=h1r-wr*h2r+wi*h2i
       data(i4)=-h1i+wr*h2i+wi*h2r
       wtemp=wr
       wr=wr*wpr-wi*wpi+wr
       wi=wi*wpr+wtemp*wpi+wi
    END DO
    if (doInverse) then
       h1r=data(1)
       data(1)=c1*(h1r+data(2))
       data(2)=c1*(h1r-data(2))
       call DFT(data,.TRUE.)
    else
       h1r=data(1)
       data(1)=h1r+data(2)
       data(2)=h1r-data(2)
    endif
  END SUBROUTINE RealFT
!C  (C) Copr. 1986-92 Numerical Recipes Software D041&0(9p#3.

  SUBROUTINE DFT(data,flagInverse)
    USE Parameters
    REAL(d),INTENT(INOUT) :: data(:)
    LOGICAL,OPTIONAL :: flagInverse ! .TRUE. for inverse DFT
    LOGICAL :: doInverse=.FALSE.
    INTEGER :: i,istep,j,m,mmax,n
    INTEGER :: isign=1
    REAL(d) :: tempi,tempr
    REAL(d) :: theta,wi,wpi,wpr,wr,wtemp
    n=SIZE(data)
    IF(PRESENT(flagInverse)) doInverse=flagInverse
    IF (doInverse) isign=-1
    j=1
    DO i=1,n,2
       if(j.gt.i)then
          tempr=data(j)
          tempi=data(j+1)
          data(j)=data(i)
          data(j+1)=data(i+1)
          data(i)=tempr
          data(i+1)=tempi
       endif
       m=n/2
       DO WHILE((m.ge.2).and.(j.gt.m))
          j=j-m
          m=m/2
       END DO
       j=j+m
    END DO
    mmax=2
    DO WHILE(n.gt.mmax)
       istep=2*mmax
       theta=2.0_d*PI/REAL(isign*mmax,d)
       wpr=-2.0_d*sin(0.5_d*theta)**2
       wpi=sin(theta)
       wr=1.0_d
       wi=0.0_d
       DO m=1,mmax,2
          DO i=m,n,istep
             j=i+mmax
!DEREK             tempr=sngl(wr)*data(j)-sngl(wi)*data(j+1)
!             tempi=sngl(wr)*data(j+1)+sngl(wi)*data(j)
             tempr=wr*data(j)-wi*data(j+1)
             tempi=wr*data(j+1)+wi*data(j)
             data(j)=data(i)-tempr
             data(j+1)=data(i+1)-tempi
             data(i)=data(i)+tempr
             data(i+1)=data(i+1)+tempi
          END DO
          wtemp=wr
          wr=wr*wpr-wi*wpi+wr
          wi=wi*wpr+wtemp*wpi+wi
       END DO
       mmax=istep
    END DO
  END SUBROUTINE DFT
!C  (C) Copr. 1986-92 Numerical Recipes Software D041&0(9p#3.

    END MODULE FFT
